A mesoscopic approach for multi-phase flows in nano-corrugated channels 
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An approach based on a lattice version of the Boltzmann kinetic equation for describing multi- 
phase flows in nano- and micro-corrugated devices is proposed. We specialize it to describe the wet- 
ting/dewetting transition of fluids in presence of nanoscopic grooves etched on the boundaries. This 
approach permits to retain the essential supra-molecular details of fluid-solid interactions without 
surrendering -actually boosting- the computational efficiency of continuum methods. The meso- 
scopic method is first validated quantitatively against Molecular Dynamics (MD) results of Cottin- 
Bizonne et al. [Nature Mater. 2 237 (2003)] and then applied to more complex situations which are 
hardly accessible to MD simulations. The resulting analysis confirms that surface roughness and 
capillary effects may conspire to promote a counter-intuitive but significant reduction of the flow 
drag with substantial enhancement in the mass flow rates and slip-lengths in the micrometric range 
for highly hydrophobic surfaces. 
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The motion of fluids at the micro and nanoscalc is con- 
trolled by the competition of dissipative effects and pres- 
sure drive. The weakness of inertia in the microworld im- 
plies that it is increasingly difficult to push fluids across 
micro/nanoconfined geometries, as their surface/volume 
ratio is made larger and larger. An obvious consequence 
is that the dynamics of microfiows is crucially affected 
by the interaction of the fluid with the confining solid 
boundaries. Information on these interactions is usually 
conveyed into the formulation of proper boundary con- 
ditions for the fluid flow. In particular, slippage prop- 
erties (see 0] for a recet review) have been reported in 
experiments and in molecular dynamics simulations, de- 
pending on the thermodynamical and wetting properties 
of the boundary (contact angle) and on the surface ge- 
ometry [J H |1 II H Q, HEI. A fundamental ques- 
tion arises as to whether the fluid really slips over the 
surface, or rather the indirect measurements based on 
pressure/mass flow rate relations or surface force appa- 
ratus reflect an apparent slip arising from surface inho- 
mogencitics or complex interface with additional physics. 
Indeed, it has been argued that a gas layer at the inter- 
face would alter the fluid dynamics in the bulk, leading 
to a mass flow rate increase even in the presence of pure 
no-slip [3 [H [HIH HI!. This hypothesis is supported 
by the observation of nanobubbles trapped on the surface 
|l3j and by a decreasing apparent slip length as the fluid 
is degased Q. 

The aim of this article is to discuss the complex effects, 
at the hydrodynamical scales, induced by the surface 
wetting properties in presence of complex geometries in 
micro- and nano-devices. The results can be summarized 
in two main points. First, we provide neat evidence that 



the physics of the boundary conditions is quantitatively 
reproduced by modeling the fluid at mesoscopic level, by 
means of a minimal version of the Boltzmann equation, 
i.e. the Lattice Boltzmann Equation (LBE) [lallq . This 
result is obtained by performing a quantitative compari- 
son of a "finite-volume" dewetting transition a gain st re- 
cent Molecular Dynamics simulations (MD) Q, h7l. Far 
from being a mere technicality, this result opens the way 
to numerical investigations at spatial and time scales 
much larger than those currently available in most MD 
simulations. Second, we extend the MD results by inves- 
tigating the critical dependency of the mass flow rate on 
the degree of roughness at constant bulk pressure. The 
simplest LBE reads as follows [l9j : 



f i (x + c i At,t + At) -fi(x,t) = 
-u>At[fi(x,t) - f( eq \x,t)}+ FiAt, 



(1) 



where fi(x,t) is the probability of finding a particle 
at site x at time t, moving along one of the b-th lat- 
tice direction defined by the discrete speed C; with i = 
1, and At is the time unit. The left-hand side 

of stands for molecular free-streaming, whereas the 
right-hand side represents molecular collisions. These 
are expressed through a simple relaxation towards local 
Maxwcllian equilibrium fj-^ in a time lapse of the order 



of t = 



Finally the term Fi represents a volumet- 



ric body-force, which can be tailored to produce highly 
non-trivial macroscopic effects, such as phase-transitions. 
Non-ideal effects, leading to two-phase flows, are modeled 
through a self-consistent force term: 



F(x, t) =GbY2 w i^( x i t)ip(x + CiAt, t)ci 



(2) 



2 



Here, ip( x ) is a phenomcnological pseudo-potential (gen- 
eralized density), ib(x, t) = ip[p(x, t)], first introduced by 
Shan and Chen |20j , Wi are normalization weights and Gb 
tunes the molecule-molecule interaction, i.e. it plays the 
role of the normalized inverse temperature, e/KT, with 
e the molecular interaction, K The Boltzmann constant 
and T the system temperature. Here we choose the stan- 
dard form ip = y/po{l — exp(— p/po)}, with the reference 
density po = 1, in lattice units. 

In spite of its simplicity, the Shan-Chen approach pro- 
vides two crucial ingredients of non-ideal fluid behavior: 
a non-ideal equation of state and a non-zero liquid-vapor 
surface tension, o\ v . Both features are encoded in the 
expression of the non-ideal momentum flux tensor Pjy. 
In the hydrodynamic limit, the LBE equations Q-© 
can be shown to evolve according to the Navier-Stokes 
equations 0] with a pressure tensor P : 



Pkj = 



1„4 



-^CgQbdkipdjip. 
The equation of state in the bulk is P 



c 2 s P 



hi 
(3) 



where we recognize a non-ideal contribution on top of 
the ideal equation of state P = pc 2 with c 2 the sound 
speed velocity. This equation of state supports a phase- 
transition at a critical density p c /po = hi(2), whenever 
the coupling strength exceeds (in magnitude) the critical 
value Qc. = —4.0. Concerning the velocity field, follow- 
ing |l8(, we set it to be zero at the boundary by us- 
ing bounce-back boundary conditions while, for the case 
when slippage properties are apriori imposed, different 
boundary conditions should be used H^,|2^|. Next, the 
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FIG. 1: Left: homogeneous roughness. A groove with depth 
h = 33Aa: and width H = L x — a (with a — lOAx) is intro- 
duced on the bottom wall and periodic boundary conditions 
are assumed along x. In this configuration, the presence of 
vapor pockets inside the groove changes the "effective" bound- 
ary conditions felt by the bulk fluid, with a net decrease of 
drag when a pressure drop is applied. Right: channel with 
heterogeneous roughness. Two grooves of width H\ = 40Aa; 
and H2 = 70 Ax are present. The two grooves are filled sep- 
arately at different values in the pressure/density diagram. 
The lattice spacing corresponds to Ax ~ 0.3nm. 



density field p{x) is assumed to match a given value 
ip w = ip(Pw), where p w should be regarded as a free- 
parameter related to the strength of the fluid/solid in- 



teractions. In 0| we have shown that by imposing 
the condition of mechanical equilibrium, djPkj = 0, of 
the contact line separating the liquid, vapor and solid 
phases, one can compute analytically the contact angle, 
6{Gb,Pw)- Notwithstanding their inherently mesoscopic 
character, the parameters Gb and p w carry no less phys- 
ical content than their atomistic counterparts (relative 
strength of the attractive to repulsive interactions in MD 
simulations 0). In order to study the wetting/dewetting 
transition on micro/nano-patterned surfaces (see Fig. ^) 
we have integrated the LBE equation JTJ in a 2D Lattice 
using the nine-speed 2DQ9 model (b = 9), one of the 
most used 2D-LBE scheme, due to its superior stability 
[2ll Eif . We have used the two geometries described in 
Fig. ^ All simulations have been performed at fixed "in- 
verse temperature" Gb = —6.0 i.e. beyond the critical 
value Gb = —4.0. The relaxation parameter is taken as 
t = 0.8 in lattice units. The equation of state delivers 
the corresponding values of the liquid and gas density, 
pi = 2.65 and p g = 0.07 respectively, while the surface 
tension is ai v = 0.105 ± 0.002 in lattice units. The value 
of the lattice spacing in physical units is obtained by 
matching the physical value of the liquid-vapor surface 
tension, with the one measured on the lattice, via the 
dimensional relation: u\^ ys — KT / (Ax) 2 ai v . For wa- 
ter/vapor at T — 540°, where the density ratio is close to 



the values of our simulation, wc have a 



phys 



0.022iV/m 



which yclds Ax ~ 0.3nm. This value is comparable with 
the atomic range of the Lennard- Jones potential used in 
MD. Molecular dynamics simulations have recently re- 
ported that the concerted effects of wetting phenomena 
and nano-corrugations can lead to a fairly substantial 
reduction of mechanical drag Q, . Specifically, the au- 
thors in [H.ll7j| consider a nanometric channel flow with a 
regular sequence of longitudinal or transverse steps (with 
respect to the mean flow) along the solid wall of the chan- 
nel, and show that, under suitable thermodynamic and 
geometric conditions, the presence of the steps triggers 
the formation of a gas film in the grooves within the ob- 
stacles. The liquid can then slide-away over the gas film, 
thereby experiencing a significantly reduced mechanical 
drag. Such phenomenon may occur only at a critical 
pressure drop between liquid and vapor phase, of the or- 
der of the capillary pressure, P cap , given by the estimate: 
0: 



P, 



2ai v cos(9) 



cap 



(4) 



L x — a 

In Fig. |3 we validate the LBE by a direct comparison 
with the MD results published in Q, 0] with the same 
geometry and with comparable contact angle. In partic- 
ular, we show the Pressure drop between the bulk liquid 
and vapor phases APi v = Pi — P v , at changing the dis- 
tance d (see Fig. ^) at fixed total mass, i.e. at changing 
the average density. As one can clearly see, the agreement 
between LBE and MD is quantitative. The plateaux ob- 
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FIG. 2: Validation of our LBE simulation against MD results 
Ht11 . The LBE is done with the same aspect-ratios of MD. The 
dimensionless normalized pressure drop, APi v h/ai v between 
the two bulk phases, is shown as a function of the normal- 
ized distance d/L x (see Fig. 0. LBE (□) results have been 
obtained with a contact angle 6 = 160°, which is consistent 
with the contact angle measured in MD (o) simulation (see 
Fig. 2 of lrH|L The two insets represent the density configura- 
tion at the onset of the wetting/dewetting transition (right) 
and for a wetted configuration (left). The plateaux in the 
pressure curve defines the capillary pressure, P cap - 



served for AP;„ in the range 0.9 < d/L x < 1.05 corre- 
sponds to the pressure/density values at which the fluid is 
invading the corrugation, forming an interface which does 
not yet touch the bottom of the groove. This corresponds 
to the capillary pressure, P cap - Reducing further d, i.e. 
increasing the average density, a change of concavity for 
0.8 < d/L x < 0.9, is observed. This range corresponds 
to values such that the interface starts to touch the bot- 
tom of the groove, adjusting its pressure/density in such 
a way as to minimize the free energy in the presence of 
the new liquid-vapor-solid interface. The agreement of 
LBE with the capillary pressure J1J is checked in detail 
in Fig. |31 where we report the change of the pressure 
diagram with a changing distance d, for three different 
corrugation values, L x — a. In the inset of Fig. |21 we 
extract the value of a\ v cos(#) from the the slope of the 
observed plateaux vs Ax/(L X — a). The value of cos(6*) 
is then obtained by estimating the surface tension, o~i v , 
through Laplace's law for a droplet in equilibrium with 
its saturated vapor. The agreement of the contact angle 
measured in this way, 9 = 158° ± 6°, with the analytical 
estimate, = 160°, obtained in 0] by imposing the me- 
chanical equilibrium condition of the contact line is very 
satisfactory. 

The comparison shown in Fig. [2]and Fig. [^demonstrates 
the first result discussed in this paper, namely that the 
model introduced in |l8l | captures the correct interplay 
between roughness and wetting effects. 
Even more complex behavior is observed for heteroge- 
neous nano-corrugations, with the simplest case shown 



FIG. 3: Pressure variations at changing d for various rough- 
nesses, L x = 50 Ax (A), L x = 60Aa; (□), L x = 80Ax (o) with 
a — 10 Aa;. In the inset, we show the relation The slope 
is given by 2a iv cos(0) = -0.196 ±0.006 (LB units). With our 
surface tension ai v = 0.105 ± 0.002 (LB units), this implies a 
best estimate of = 158° ± 6°. 



in the right panel of Fig. ^ In this case, one has 
two characteristic groove sizes, Hi and H^, and hence 
two corresponding critical capillary pressures. The pres- 
sure/density diagram for this heterogenous corrugation is 
shown in Fig. 0] where the two plateaux corresponding 
to the two capillary pressures coexist. This device may 
be considered a "smart" two-state surface, whose wetting 
properties and mass throughput (under the application 
of a pressure gradient ) may be tuned by changing the 
bulk pressure. 

The dynamical response of the micro-channel is investi- 
gated by applying a constant pressure gradient. In Fig.0 
we show the mass flow rate as a function of the degree of 
roughness, a/(L x — a), for a fixed groove depth h and at a 
given normalized pressure drop, APi v h/ai v = 0.75. The 
main result here is the presence of a transition at a criti- 
cal roughness, a/(L x — a), where the mass flow rate starts 
to increase with respect to the perfect wetting situations 
reaching as much as 100% gain. The strong dynamical 
effect of the gas-layer can be quantified in more detail by 
inspecting the velocity and momentum profiles along the 
vertical direction without the vapor layer (fully-wetted 
configuration) and with a thin vapor layer starting to 
accumulate close to the bottom of the groove. This is 
shown in the inset of Fig. [S] where the local slip length 
is also depicted by extrapolation of the bulk profile in- 
side the wall. As soon as a vapor layer is formed, the 
local slip-length ramps-up, reaching values of the order 
of the channel height 45Aa; ~ 15ran. Even larger values 
can be measured close to the dewetting transition, where 
a well-developed vapor layer is formed inside the groove 
(the so-called super-hydrophobic regime). 

Summarizing, we have shown that an extension of lat- 
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FIG. 4: Pressure variations at changing d, with fixed total 
mass for the case of heterogeneous roughness (see right panel 
of Fig.0. The two plateaux correspond to the case in which 
the liquid is starting to invade the widest groove (right) and 
when it is invading also the thinnest (left). 



tice Boltzmann equation for non-ideal fluids, can quan- 
titatively account for the concerted effects between wet- 
ting phenomena and geometrical irregularities. In par- 
ticular, the presence of nano/micro-irrcgularitics in the 
flow geometries leads to sizeable effects with respect to 
the infinite volume liquid-gas transitions, as well as to a 
significant reduction of mechanical drag on the flowing 
fluid. The consequence of our results is twofold: from a 
theoretical perspective, it indicates that drag-reduction 
via geometry-induced wetting transitions is a non-specific 
phenomenon. On the practical side, the present LBE ap- 
proach offers the opportunity to perform very efficient 
numerical simulations of complex micro/nanofluidic phe- 
nomena at scales of direct experimental relevance, which 
are hardly accessible to atomistic simulations. This will 
open the way to the systematic optimization of microflu- 
idic devices via computer simulation. Useful discussions 
with J.-L. Barrat, X. Shan and S. Troian are kindly ac- 
knowledged. 



FIG. 5: Mass flow rate normalized to fully wetted case vs 
the effective roughness a/{L x — a). The bulk pressure is 
fixed to be AP;„ = 0.75ai v /h. A critical roughness (verti- 
cal line) is given by the estimate of the capillary pressure, 
0.75 7 
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